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Abstract 

We study the qualitative properties of population cycles in a predator-prey system where genetic vari- 
ability allows contemporary rapid evolution of the prey. Previous numerical studies have found that prey 
evolution in response to changing predation risk can have major quantitative and qualitative effects on 
predator-prey cycles, including: (i) large increases in cycle period, (ii) changes in phase relations (so that 
predator and prey are cycling exactly out of phase, rather than the classical quarter-period phase lag), and 
(iii) "cryptic" cycles in which total prey density remains nearly constant while predator density and prey 
traits cycle. Here we focus on a chemostat model motivated by our experimental system 1 14 41 ] with algae 
(prey) and rotifers (predators), in which the prey exhibit rapid evolution in their level of defense against 
predation. We show that the effects of rapid prey evolution are robust and general, and furthermore that 
they occur in a specific but biologically relevant region of parameter space: when traits that greatly reduce 
predation risk are relatively cheap (in terms of reductions in other fitness components), when there is co- 
existence between the two prey types and the predator, and when the interaction between predators and 
undefended prey alone would produce cycles. Because defense has been shown to be inexpensive, even 
cost-free, in a number of systems I4lll6ll42l . our discoveries may well be reproduced in other model sys- 
tems, and in nature. Finally, some of our key results are extended to a general model in which functional 
forms for the predation rate and prey birth rate are not specified. 
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1 Introduction 



The potential for rapid evolutionary change is well documented: in some natural and experimental settings, 
trait evolution and population dynamics are observed to occur on similar time scales [20 22 29 38 1. Though 
laboratory demonstrations of the interactions between trait and population dynamics still remain rare, there 
are a few direct examples 1131 l27l and also examples where these interactions are inferred via modeling 
(e.g., 1110114111431 ). Furthermore, documented cases where the effects of rapid evolution are observed in a 
changing natural environment are increasing (e.g., ll5l 1121 1171 1191 121 1 1281 1321 . and for reviews on the topic, 
j6]|20||33| 44 1). Among the fishes, rapid evolution has been observed or inferred in natural populations 
of atlantic cod (Gadus morhua), which evolved towards reproductive maturity at earlier ages and smaller 
sizes under heavy size-dependent selection (fishing) H28I . Within four hatchery generations, Chinook salmon 
{Oncorhynchus tshawytscha) raised in hatcheries to supplement natural populations evolve life-history traits 
(smaller eggs) that reduce odds of survival in the wild |21 1. 

Where prey are under strong selection to avoid predation (because the risk of getting eaten is very strong nat- 
ural selection), prey genetic diversity yields the capacity for rapid evolution of resistance to predation, akin 
to the rapid evolution of microbial pathogens in response to antibiotics. However, prey defensive strategies, 
be they physiological, behavioral, or developmental, may themselves exact demographic costs from the prey 
equivalent to or greater than that of direct consumption by predators 1 3 1 1 . The energetic costs of defense traits 
in "trait mediated interactions", for example the development of armored spines in Daphnia when exposed 
to chemicals from fish, reduce lifetime fitness by routing energy to defense which would otherwise go into 
progeny Q. Thus, by focusing exclusively on changes in population numbers without considering changes 
in the properties (and associated costs) of individuals comprising the populations, conventional models of 
population and community dynamics may give us, at best, only half the story in a system where prey adapt 
to predation risk. The capacity for rapid evolution may expand the conditions permitting coexistence of 
predator and prey. For example, within the ecosystem of the human body, rapid prey evolution may allow 
viruses (prey) to escape detection by immune cells (predators) by altering a few surface proteins, or allow a 
pathogenic strain of bacteria to evolve so that it is no longer sensitive to a widely-used antibiotic. Hetero- 
geneity in prey edibility may also shift prey population regulation from top-down to bottom-up by providing 
a haven from predation in the form of the less edible prey. 

Our experimental system is a predator-prey microcosm with rotifers, Brachionus calyciflorus, and their algal 
prey, Chlorella vulgaris, cultured together in nitrogen-limited, continuous flow-through chemostats. In prior 
studies, we have shown that coexistence of edible and inedible prey types (genetic variation in the algal prey) 
allows the prey to evolve in response to temporally variable selection due to predation pressure from the rotifer 
predator, and nutrient limitation at high prey densities. While our earlier studies did not specifically track 
changes at the genotypic level, recent work on our system explicitly identifies two competing algal strains, 
and tracks changes in their densities as the algal population evolves under predation pressure 1 27 1 . Evolution 
in the prey can lead to "evolutionary" cycling 1341 1411 . where the predator and prey exhibit extended, out- 
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of-phase population cycles (Figure^V), or in some circumstances, the odd phenomenon of "cryptic cycles", 
where the predator alone exhibits regular population cycles but the prey appear to remain in steady state 
(Figure In cryptic cycling, densities of edible and inedible prey cycle out of phase with each other, 
driven by changes in predator abundance, in such a way that total prey density remains nearly constant 1 4-3 1. 
These dynamics are not unique to the organisms in this system: we have observed evolutionary cycles in 
a chemostat system comprised of rotifers cultured with the flagellated algae Chlamydomonus, and cryptic 
dynamics have been observed in bacteria-phage microcosms II 101 1431 . We are motivated here by just these 
sorts of perplexing experimental results from our system, and by the close match between our experimental 
data and model simulations. 

Understanding the potential effects of rapid evolution on the dynamics of natural ecosystems is critical to 
predicting how populations will adapt to a changing environment. Populations in the wild today face un- 
precedented stress from habitat loss or degradation, harvesting pressure, species introductions and climate 
change. In addition, otherwise well-intentioned attempts at conservation or management often fail to take 
into account the potential for rapid Darwinian responses to intervention |24|. Thus before conclusions based 
on laboratory systems or manipulated natural systems are applied to the natural world, we must ask if the 
conclusions are likely to be robust: are they limited to the special conditions present in the experimental 
systems, or should we expect to see them in a broad range of conditions in nature? The present contribution 
is an attempt, using theory, to answer the questions: how general is the phenomenon of evolutionary cycling 
in predator-prey systems, under what circumstances might these dynamics be observed, and what are the 
implications of this type of phenomenon for natural systems? 

2 The model 

Our model is based on an experimental predator-prey microcosm with rotifers, Brachionus calyciflorus, and 
their algal prey, Chlorella vulgaris, cultured together in a nitrogen-limited, continuous flow-through chemo- 
stat system. This system was first described by Fussmann et al. (l4i|, further characterized by Schertzer et 
al. [34 1 and Yoshida et al. 141 1 1421 . and equilibrium properties studied by Jones and Ellner |23|. Brachionus 
in the wild are facultatively sexual, but because sexually produced eggs wash out of the chemostat before 
offspring hatch, our rotifer cultures have evolved to be entirely parthenogenic 1 15 1. The algae also reproduce 
asexually |30|, so evolutionary change in the prey occurs as a result of changes in the relative frequency of 
different algal clones. 

We use a system of ordinary differential equations to describe the population and prey evolutionary dynamics 
in the experimental microcosms H4 1 1 1231 . Genetic variability and thus the possibility of evolution in the prey 
is introduced by explicitly representing the prey population as a finite set of asexually reproducing clones. 
Each clone is characterized by its palatability p, which represents the conditional probability that an algal 
cell is digested rather than being ejected alive, once it has been ingested by a predator 1 27 1 . 
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A: Evolutionary Cycling 




Figure 1 : A, "Evolutionary" cycling. Chlorella are shown in solid line, and the rotifer predator is shown dashed. Phase 
relations are nearly out-of-phase, unlike the quarter-phase shift seen in classic predator-prey cycles. B, Cryptic Cycling. 
Initially the system exhibits classic predator-prey cycles, which would be expected when a single (edible) prey type is 
dominant. At about day 50 the system switches to cryptic cycles, which would be expected if a highly defended (inedible) 
type with low cost for defense arose by mutation. A switch from classic to cryptic cycles when a defended type arises by 
mutation has been documented in bacteria-phage chemostats 1 10 1. 



The model consists of three equations for the limiting nutrient and rotifers, plus q equations for q prey clones. 
In the following equations, N is nitrogen (p, mol per liter), C, represents concentration of the i th algal clone 
(10 9 cells per liter), where i = l,2,....,q. Here we limit the number of clones to two, for reasons discussed 
below. R and B are the fertile and total population densities, respectively, for the predator Brachionus (indi- 
viduals per liter). Fertile rotifers senesce and stop breeding at rate X; all rotifers are subject to fixed mortality 
m. The parameters XciXb 816 conversions between consumption and recruitment rates (additional model 
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parameters are defined in Table 1). 



dN q 

— = 8{N I -N)-^F c>i {N)C i 

i=l 

-T = XcFcAN)Q - Pi F b {Ci)B - SC, 
dt 

d_R * CD 

dt 



=1 

^=Xbt,PiF b (Q)R-(8 + m)B 

where 

F cJ (N)=p c N/(K c ( Pi )+N) 

and 

2 

F b {C i )=GC i /(K b + Y J C i pi) (2) 

/=i 

are functional response equations describing algal and rotifer consumption rates, respectively, and where 
p c = (OcPc/ e c . Equation is derived from the predator's clearance rate G (the volume of water per unit time 
that an individual filters to obtain food). We assume that clearance rate is a function of the total prey food 
value: 

G= % . (3) 

i=l 

That is, lower prey palatability results in the predators increasing their clearance rate, exactly as if prey were 
less nutritious. We also considered a model in which clearance rate depends only on the total prey density, 
but it could not be fitted as well to our experimental data on population cycles. Elsewhere 1411 1271 we have 
used a more complicated expression for G in order to fit experimental data more accurately, but using 
does not change the model's qualitative behavior. 

The cost for defense against predation is a reduced ability to compete for scarce nutrients 11421 1231 1771 . We 
model this by specifying a tradeoff curve 

K c (p)=K c + a 2 (l-p) a i. (4) 

Here K c > is the minimum value of the half-saturation constant, a\ > determines whether the tradeoff 
curve is concave up versus down, and a-i > is the cost for becoming completely inedible (p = 0). 
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Table 1 : Parameter estimates for the Chlorella-Brachionus microcosm system. Set: adjustable parameters 
set by the experimenter. TY: Unpublished experimental data (Yoshida et al, in prep.) Fitted: Estimated 
by numerically optimizing the goodness-of-fit between model output and data on total prey and predator 
population dynamics from two experiments (originally reported by Fussmann 1 14|) in which regular cycles 
occurred. 



Parameter 


Description 


value 


Reference 


Nj 


Limiting nutrient cone, (suplied medium) 


80/1 mol N/7 


Set 


8 


Chemostat dilution rate 


variable (d) 


Set 


V 


Chemostat volume 


0.33Z 


Set 


Xc 


Algal conversion efficiency (10 9 cells//imol N) 


0.05 


|14| 


Xb 


Rotifer conversion efficiency 


ps 54000 rotifers/10 9 algal cells 


Fitted 


m 


Rotifer mortality 


0.055/d 


|14| 


I 


Rotifer senescence rate 


0.4/d 


|14| 


K c 


Minimum algal half-saturation 


4.3 ii mol N // 


m 


K b 


Rotifer half-saturation 


0.835 x 10 9 algal cells // 


TY 


A 


Maximum algal recruitment rate 


3.3/d 


TY 




N content in 10 9 algal cells 


20.0 Li mol 


|14| 


£c 


Algal assimilation efficiency 


1 


|14| 


G 


Rotifer maximum consumption rate 


5.0 x 10- 5 l/d 


TY 


CC\ 


Shape parameter in algal tradeoff 


variable, (X\ > 


Fitted 


«2 


Scale parameter in algal tradeoff 


variable, cfe > 


Fitted 



3 Characteristics of the model under simulation 



A system of q > 2 prey types invariably collapses to at most two types in the presence of a predator: either 
a single clone that outcompetes all others, or a pair of very different clones (one very well defended and the 
other highly competitive) that together drive all intermediate prey types to extinction |41 23 1. Only the latter 
case is of interest here, because with a single prey type there is no prey evolution. We thus consider here a 
system of two extreme prey types in the presence of a predator. 

Two system parameters can be experimentally varied: the dilution rate 8 (fraction of the culture medium that 
is replaced daily) and the concentration of the limiting nutrient in the inflowing medium, N]. Fussmann et 
al. [14 1 showed that 5 is a bifurcation parameter: in both the real system and the model, the system goes 
to equilibrium at low dilution rates, limit cycles at intermediate dilution rates, and again to equilibrium at 
high dilution rates. Further increases in 8 lead to extinction of the predator. Toth and Kot |39| proved that 
the same bifurcation sequence occurs in chemostat models with an age-structured consumer feeding on an 
abiotic resource (for our experimental system, this would be rotifers feeding on externally-supplied algae that 
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could not reproduce within the chemostat). 

The prey vulnerability parameter p is also a bifurcation parameter. In the following discussion, we define 
evolutionary cycles as both prey types coexisting and exhibiting long-period cycles (period 2CM10 days), 
with the predator and total prey abundance almost exactly out-of -phase with each other. Predator-prey cycles 
are shorter (6-12 days), display the classic quarter-period phase offset between predator and prey, and involve 
one prey type cycling with the predator. In addition, both prey may survive and coexist with the predator at 
an evolutionary equilibrium, or one prey type may be driven to extinction while the other goes to equilibrium 
with the predator. 

Single prey model Figure|2j\ shows the dynamics of the single prey model as a function of prey palatability 
p and dilution rate 8. Parameters giving single -prey predator-prey cycles are indicated by open circles, 
and elsewhere the system goes to equilibrium. At low p values (up to 0.4-0.6, depending on the predator 
conversion efficiency %b) the system goes to equilibrium at all dilution rates. As p increases there is a 
bifurcation and short, low amplitude predator-prey cycles are observed, initially for the narrow range of 
dilution rates. When p is higher, oscillations grow in amplitude and increase very slightly in period, and 
cycling occurs over a larger range of dilution rates. The cycles always exhibit classic predator-prey phase 
relations. 

Two prey models Figure |2j3 shows dynamics of the two prey model as a function of the dilution rate 8 
and the trait value p of the defended prey type (the model is scaled so that the undefended type has p = 1). 
Using the parameter values listed in Table 1, extended evolutionary cycles (closed circles) initially appear for 
all dilution rates (0.2 < 8 < 1.3 at pi very small {pi 0.01). As p\ increases, evolutionary cycling occurs 
for a diminishing range of dilution rates. By pi = 0.2, cycling vanishes and instead the defended prey is 
in equilibrium (Figure |2j3 , white space) or the two prey types are in an evolutionary equilibrium with the 
predator (Figure[2j3 , crosshatching). As pi increases further (0.4 <p\ < 0.6), there is another bifurcation and 
the system, comprised of the defended type and the predator, begins to exhibit predator-prey cycles (Figure 
|2jj, open circles). From this point on the system behaves as if it were dominated by the defended type (see 
above), until pi has increased to the point that the two prey types are almost identical. At that point there are 
predator-prey cycles with both prey types present (closed circles) but these appear to be very long transients 
rather than indefinite coexistence: one or the other prey type, depending on the dilution rate, is slowly driven 
to extinction by its competitor. 

Eliminating predator age structure Panels C and D in Figure[2]shows model dynamics without age struc- 
ture in the predator. Age structure is removed by setting A = in Q, with all other parameters unchanged. 
As seen in Figure |2jl!, the single prey model without age structure exhibits dynamics very similar to those 
in|2}, where age structure is included. Predator age structure is generally stabilizing in this model because 
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senescent rotifers are a resource sink, eating prey without converting them to offspring. This effect is most 
pronounced at low values of 8 because senscent rotifers then spend more time in the chemostat before getting 
washed out. Omitting age structure is therefore destabilizing: it permits cycles with better defended prey 
(lower p) and eliminates entirely the stability at very low 8 for nearly all p values. Similarly, simulations 
of the two-prey model show that eliminating predator age structure shifts the region of (p,8) values giving 
evolutionary cycles to higher dilution rates, and eliminates the stabilization at very low 8, but otherwise the 
bifurcation diagram is unchanged. Given these similarities in model behavior, we may simplify model by 
eliminating predator age-structure without changing the properties of interest for this paper. 



4 Rescaling the model 

We now simplify the model (11 12t by rescaling and further reducing its order. Based on the simulation results 
above we omit age structure in the predator, which is now represented by the one state variable B. We also 
assume that predator mortality m is negligibly small relative to washout at the dilution rates 8 of interest, and 
set m = 0. The model then becomes: 

dN Jf/A7 » t\ £ NQ 



dt 

dCi 
dt 
dB 
dt 



= C, 
= B 



XcP, 

Xb 



N 



G Pi B 



K c ( Pi )+N (K B +ZPiCi 



-8 



(5) 



-5 



We order the prey types so that p \ and pi correspond to the defended and vulnerable prey types, respectively, 
P\ < Pi- The cost for defense is reduced ability to compete for scarce nutrients, so K c {p\) > K c (p2). 



To rescale the model we make the following transformations: 

B 



S=- Q 



y - 



XbG 
8 ' 



XcPc 



k h = 



K B 



XcXbNi 1 ° 8 ' 8 ' - XcNi 

The half-saturation constants for each of the two prey types are transformed as follows, 



t = 8t 



k x =K c { Pl )/N h k 2 =K c {p 2 )/Ni. 



(6) 



(7) 



Substituting these into gives: 



y = y 



-s-s] 

i 


^k + S 


mS 




gpty 


_ki + S 


(kb + Q) 


\ gQ 


-1 




Ub + Q 





- 1 



(8) 



where 



Q=P\X\ +p 2 X2- 



(9) 
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Table 2: Estimates of rescaled model parameters for the Chlorella-Brachionus microcosm system. 



Parameter 


Description 


Value 


m 


Algal maximum per-capita population growth rate 


3.3/5 


k\M 


Algal half-saturation constants for nutrient uptake 


0.054 


8 


Predator maximum grazing rate 


2.55/5 


h 


Predator half-saturation constant for prey capture 


0.21 



Table 2 gives values of the rescaled model parameters corresponding to the parameter estimates in Table 1 . 

We can reduce the dimension of the system further by letting TL = S-\-x\ +X2+y. From (|8} we have E = 1 — E, 
so E(f ) — > 1 quickly as t — > °°. Thus, to study the long-term dynamics of 0, we may consider the dynamics 
on the invariant set E = 1 . Then S{t) = 1 —x\{t) — X2 (f ) —y(t), and (|8j reduces to 

m(l -X-y) 



X{ X; 



y = y 



gpiy 

ki + (l-X-y) (k b + Q) 

gQ 



i= 1,2 



k b + Q 



- 1 



(10) 



where X = x\ +X2- 



5 Analysis 

Our goals in this section are to find the conditions under which two prey types can coexist, to determine when 
coexistence is steady-state versus oscillatory, and to characterize the period of cycles and the phase relations 
during oscillatory coexistence and during transients when one type is decreasing to extinction. Throughout 
this section we consider the reduced model JlOi . For local stability analysis it is useful to note that the model 
has the form 

Xi = x i n(xi 1 x 2l x i ) i = 1,2,3 (11) 

with X3 = y. It follows that at any equilibrium where the x, are all positive (and hence the r,- are all 0) the 
Jacobian matrix J has entries 

J{iJ)=x iJ p- (12) 

dXj 

with the tilde indicating evaluation at the equilibrium with all x, present. It is also useful for local stability 
analysis that the determinant of d!2i is always negative unless p\=p2 (AnnendixlDl. 
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5.1 Dynamics of a one-prey system 



We need first some properties of the one-prey model 

w(l — x — y) 



gpy 



y = y 



k+(l— x— y) (kb+px) 



1 



(13) 



_k b + px 

This is a standard predator-prey chemostat model and its behavior is well-known, so we summarize here only 
the results that we will need later; see e.g. [36 1 for derivations and details. 

In the absence of predators, the steady state for this system is Eq = (1 — A,0), where 

k 



A = 



m — 1 



(14) 



A is the scaled concentration of limiting nutrient at which prey growth balances washout rate, so that x = 0. 
Similarly, steady state densities for each prey type in a predator-free two clone system are 



Ei={\— A;,0) where A,- = 

The steady state for the prey in the presence of the predator is 

h 



m— 1 



(15) 



x c is the prey density at which the predator birth and death rates are equal. The model Jl 3I > has an interior 
equilibrium point E c = (x c ,y c ) representing predator-prey coexistence if 



A+x c < 1 



(16) 



1 36 1, and otherwise the predator cannot persist. The system then collapses to the prey by itself and converges 
to Eq. Condition ( 1 16b says that there is an interior equilibrium if the prey by themselves reach a steady state 
(1 — A) that provides enough food so that the predator birth rate exceeds the predator death rate. 



The expression for the steady state of the predator, y c , is easily obtained from (1 1 31 : 



y c 



a-x c , where a = (x c +y c ) = ^ 



y- \J {f-\mx L ) 



(17) 



with y = k+l + mx c - Similarly, the steady-state densities for the predator in a single-prey system with either 
prey type, yi, are found by substituting the steady state for the prey, x,, in place of x c and the appropriate 
half-saturation ki in place of k in d!7l >. 

We can use (I16> to derive the condition for predator-prey coexistence in terms of the prey defense trait p and 
the dilution rate 8, recalling that A and x c are both implicit functions of 8. Using J14l> and dl5> we obtain 
from (TlBl 



m(S)-l P g(8)-l 



< 1 , or p > 



l-A \g(S)-l 



(18) 
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The quantity within parenthesis above is the amount of substrate present in perfect food (undefended prey 
with p = 1). Solving Jl 8i for 8 in terms of p yields the boundary between predator extinction and stable 
coexistence in Figure [3] To the left of this line, the predator goes extinct and the equilibrium is Eq. As the 
left-hand side of the second expression in is an increasing function of p, and the right-hand side is an 
increasing function of 8, the range of p values yielding coexistence narrows as 8 increases (see Figure^}- 

As in the standard Rosenzweig-MacArthur predator-prey model, the stability condition has a graphical inter- 
pretation in terms of the nullclines. The prey nullcline is a parabola which peaks at 

. 1 



l-*+VA(m-2) 

The coexistence equilibrium is locally unstable if the peak of the prey nullcline is to the right of the predator 
nullcline (i.e., if x* > x c ). Note that a system with defended prey (p < 1) is always more stable than a system 
with fully vulnerable prey (p = 1) as reductions in p shift the predator nullcline to the right. 



From d!2t the Jacobian of Jl 3I > at E c has the form 

Jc = 



- dx 



(19) 



so E c is locally stable if the trace Tr(J c ) = *c^p is negative. Cycles emerge through a Hopf bifurcation when 
the trace becomes positive. The condition Tr(J c ) > is equivalent to the following expression for model 
(EJ: 



mk 



gp 2 yc 



> 61 JC (70) 

{k+\-x c -y c f ~ (k„+px c ) 2 ^ 

11361 . Cycles begin when the rates of change in prey substrate uptake (LHS) and in predator consumption 
(RHS) with respect to the amount of substrate present as prey (x) are exactly equal. Numerically solving (I20> 
for 8 in terms of p yields the boundary between stable coexistence and predator-prey cycles in FigurefJ] It is 
known that these cycles are stable and unique near the Hopf bifurcation, and numerical evidence uniformly 
indicates that they are always unique and attract all interior initial conditions except E c 1 36 1. 



5.2 Stability and dynamics of a two-prey system 

System dlOt has two prey types ordered so that < pi < p2 = 1 ■ We refer to prey 1 as the defended type and 
prey 2 as the vulnerable type. The cost for defense comes in the form of reduced growth rate, \/k\ < 1/&2- 

Following Abrams 1 3 1, we begin by finding the conditions for existence of an equilibrium (x\ ,X2,y) at which 
all three population densities are positive; we refer to this as a coexistence equilibrium. Setting i ll Ob to 
zero and solving gives expressions for X,Q and y in terms of model parameters (see Appendix Icl as above 
Q = p\X\ + pix% is the total prey quality, and X = x\ +X2 is the total prey density). The prey steady states 
x. i,j?2 are then 



X\ 


1 


' PiX-Q ' 


. * 2 . 


P2-P1 





(21) 
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where 

Q=\- (22) 
A coexistence equilibrium thus exists provided X > and p\X < Q < /?2^, or 

pi<f<p 2 - (23) 

Beyond the above, system dlOi does not yield tidy analytical solutions for the steady states at coexistence. 
To study how parameter variation affects coexistence, we start by graphically mapping the region where a 
coexistence equilibrium exists as a function of the defended clone's parameters, p\ and k\ (Figure|4}, without 
regard to whether or not the equilibrium is locally stable. The coexistence region also varies with 5, but 
selecting several 8 values of interest gives a general sense of how the coexistence region varies as a function 
of dilution rate. 

The lower boundary of the coexistence region occurs when the cost of defense is so high that the equilibrium 
density of the defended prey x\ drops to zero while x 2 and y remain positive. Recalling the general form il Q , 
the lower boundary is thus defined by the conditions 

n (0,x 2 ,y) = r 2 (0,x 2 ,y) = r 3 (0,x 2 ,y) = with x 2 > 0,y > 0. 

The conditions on r 2 and r 3 are solved by the steady state E 2 = (Q,x 2 ,y 2 ) for a one-prey system with only the 
vulnerable prey. The lower boundary of the coexistence region is thus defined by the condition r\ (0, x 2 ^y%) — 
0, which can be written as 

' ' ^ 1+ * 2 (24) 



k\ \-x 2 -y 2 



x 2 (m-l)-y 2 pi 



The upper boundary of the coexistence region occurs when the cost of defense is so low that the defended 
prey (at the equilibrium density) drives one of the other populations to extinction. In section [6] we show that 
for pi < p* — (5/(1 — A2), the predator goes extinct first (y — > 0) as k\ decreases, because the defended prey 
(at steady state) drives the vulnerable prey to low abundance and the defended prey is very poor food. This 
occurs atfci =k 2 (zero cost of defense). Forpi > p*, the vulnerable prey type is outcompeted by the defended 
type before k\ has reached k 2 . This boundary is therefore defined by the conditions 

ri(xi,0,y) =r 2 (xi,0,y) = r 3 (x\ ,0,y) = with x { > 0,y > 0. 

The conditions on r\ and ri are solved by the one-prey steady state E\ = (x\ ,0,yi), so the condition r 2 (f\ ,0,y~i) 
defines the upper boundary of the coexistence region for p > p*. The upper boundary of the coexistence 
region is thus the curve 

1 r 1 1 

(25) 



1 

— = mm 

h 



where <p is value of k\ that solves 



k 2 q>(xi,y~i) 
m(\~x\~y\) y\+p\x\ 



k 2 + (l-xi-yi) {pixi) 



(26) 
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noting that x\ and y~\ are functions of k\ and p\. The two segments of the upper boundary defined by j25l > 
meet at the point 



As 8 — > (with the parameter scalings in Table 2), Q J, and A 2 | 1, so p* | so there is typically an 

increasingly narrow band of p\ values for which a p\ k\ tradeoff curve lies in the coexistence equilibrium 

region (unless the tradeoff curve happens to lie exactly inside the cusp of the coexistence equilibrium region). 

As p\ — > 1, the upper and lower boundaries of the coexistence region meet at p\ = p% = l,k\ = fc 2 (Figure 
|4}. That is, if the two prey are almost equally vulnerable to predation, they can only coexist at equilibrium 
if a tiny bit of defense has a tiny cost. To prove that this occurs, we show that the point p\ = l,k\ = ki lies 
on both boundaries. At this point the two prey are identical so x\ — x 2 and y\ — y~2- The upper boundary is 
defined by r 2 (xi,0,;yi) = 0. At pi = l,k\ = k2, 

r 2 (ii,0,yi) =r 2 (x2,0,y 2 ) =^2(0,^2,^2) =0, 

thus pi = l,k\ = k2 lies on the upper boundary. The lower boundary is defined by n (O,^,)^) = ®' ^ 
p\ = l,k\ =k 2 ), 

n(0,X2,y~2) = n(0,xiyi) = r 2 (xi,0,y~i) = 0, 

which shows that p\ = l,k\ = A; 2 lies on the lower boundary. Thus both boundaries converge to k\ = k 2 as 
Pi - I- 



5.3 Local stability of coexistence equilibria 

To characterize two-prey evolutionary cycles we need to find the bifurcation curves in parameter space where 
these cycles arise. The "empirical facts" are summarized in Figure|5] based on numerical evaluations of the 
Jacobian and its eigenvalues within the coexistence equilibrium region. In Figure|5]we change the stability of 
the (predator + vulnerable prey) system by varying the value of 8, but the results are qualitatively the same if 
other parameters are varied instead (e.g., varying the prey maximum growth rates). 

The stability properties in Figure|5]explain the major qualitative features of the two-prey model's bifurcation 
diagram (Figure |2jD). To see the connection, recall that a horizontal (constant 8) slice through Fig. |2]D 
corresponds to a tradeoff curve between p\ and k\ in the panel of Figure|5]with the same value of 8. Panel 
A of Figure |3 has 8 = 1.5. When p\ is near 1, the tradeoff curve lies above the coexistence equilibrium 
region, and the defended prey type eventually outcompetes the vulnerable type. For p \ very close to 1 the 
prey types are very similar, and the vulnerable type persists for a long time. The system exhibits "classical" 
predator-prey cycles as if a single prey-type were present, even though two types are transiently coexisting. 
For p\ somewhat smaller, the vulnerable type is quickly eliminated and there are either classical cycles with 
only the defended type (open circles in Fig. |2j3), or (for lower values of p\) the defended prey type goes to 
a stable equilibrium with the predator (open triangles in Fig. |2jD). As p\ decreases further, Figure|5j\ shows 
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that the tradeoff curve enters the coexistence region in the area where the coexistence equilibrium is stable, 
so the system then exhibits stable coexistence (cross-hatching in Fig.[2f>). Finally, as p\ decreases towards 0, 
the tradeoff curve enters the area where the coexistence equilibrium is unstable, and it lies above the dash-dot 
curve marking the k\ value required for the defended prey type to invade the vulnerable prey's limit cycle 
with the predator. The system exhibits evolutionary cycles with both prey types persisting (filled circles at 
p«0inFig.|2)). 

Figure |5j\ also shows that there is a region of parameters (below the dash-dot curve) where the coexistence 
equilibrium is stable and the system therefore has coexisting attractors: a locally stable coexistence equilib- 
rium, and a locally stable limit cycle with the predator and the vulnerable prey. 

Figure|5jj, which has 8 = 1 .75, shows the same sequence of transitions as Figure|5l\, but each occurs at higher 
values of p\, reflecting the stabilizing effect of increased washout. This is reflected in Figure|2j): increasing 
8 above 1.0 shifts all the bifurcation boundaries to higher p values, but the sequence of bifurcations as p 
decreases is unchanged. However for 8 sufficiently high (panels C and D in Figure [5ji, the tradeoff curve 
lies either below the coexistence equilibrium region or within the region where the coexistence equilibrium 
is stable, so evolutionary cycles never occur. Instead, there is either stable coexistence of the two prey with 
the predator, or classical predator-prey cycles with only the vulnerable prey type. 

Evolutionary cycles are also eliminated - but for a different reason - as 8 J. in Figure|2jD. As noted above, 
as 8 i we also have p* J, 0, so unless p\ ps the tradeoff curve lies above the coexistence equilibrium region 
and only the defended prey persists with the predator, cycling at higher p\ and stable at lower p\ . Only very 
near p\ = 0, a region tiny enough to be missed by our simulation grid in Figure I3 can there be coexistence 
of both prey with the predator. 

Stability on the edges. We can gain some understanding of the patterns in Figure [5j and see that they are 
not specific to the parameter values used to draw the Figure, by examining the limiting cases that occur along 
the edges of the coexistence equilibrium region. One general conclusion (explained below) is that the bottom 
and right edges, and the right-hand portion of the top edge, all must have the same stability as the reduced 
system with the predator and only the vulnerable prey (prey type 2). However even if this system is unstable, 
there must be a region along the top edge where the coexistence equilibrium is stable. 

The Jacobian matrix that determines equilibrium stability is derived in Appendix ID1 We also show there 
that the determinant of this Jacobian is always negative at a coexistence equilibrium unless p\ = p2, so 
the coefficient cq = — det(7) in the Routh-Hurwitz stability criterion for 3-dimensional systems is always 
positive. 

Bottom and right edges: Near the bottom and right edges, the coexistence equilibrium has the same local 
stability as the (predator + vulnerable prey) subsystem (panels A and B versus C and D in Fig [5}. The 
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bottom edge is the lower limit of the coexistence equilibrium region, where x\ — > 0. The coefficients for the 
Routh-Hurwitz stability criterion (see Appendix|SJl are then 

c = -det(7) >0, ci = T 2 (J)->8 i , c 2 = -T{J)^-T 2 (27) 

where &i and T 2 are the determinant and trace, respectively, of the 2x2 Jacobian for the (predator + vulnerable 
prey) system. If this one-prey system is stable then 52 > 0, T 2 < so cq,c\ and c 2 are all positive. Moreover 
co = 0(x\) (see AppendixIDl. so whenii is small we have c\c 2 > co and the equilibrium is stable. Conversely 
if the steady state for the (predator + vulnerable prey) system is unstable, c 2 is negative so the full system is 
also unstable. 

The right edge corresponds to the cusp in the coexistence region as p\ — > 1. Near the cusp the two prey 
become increasingly similar (pi w p 2 — \,k\ ~ k 2 ). Using ill 2b . the functional forms of the r\ and the fact 
that p\ « p 2 then imply that the form of J is approximately 



Jo 



aq aq —qb 

a(l — q) a(l — q) — (l—q)b 
c c 



(28) 



where q = Xi/(x\ +JE2); even if p\ is near p 2 , it is not necessarily the case that x\ is close to x 2 . In d28t b and 
c are positive while a has the sign of dr\ /dx\ which may be positive or negative. One eigenvalue of Jq is 0, 
corresponding to the dynamics of x\ — x 2 . The others are A (a ± Va 2 — Abe) , which are also the eigenvalues of 
a single-prey system at the coexistence steady state. Thus, the two-prey system with p\ w p 2 = 1 "inherits" 
two eigenvalues from the one-prey system with p = 1. 

When the one-prey system with p=l is cyclic, the inherited eigenvalues are a complex conjugate pair. In 
the corresponding eigenvectors, the components for the two clones are identical when p\ = p 2 . This implies 
that when p\ w p 2 the eigenvector components will be similar, so the two prey types cycle almost exactly in 
phase. The period of these oscillations is determined by the inherited eigenvalues, so it is close to the period 
of the corresponding one-prey system. 

When the one-prey system is stable, the Routh-Hurwitz criterion (Appendix |SJ|, using Jq to approximate 
trace(7) and T 2 (J) and the fact that det(/) < for p\i^ p 2 , implies that the full system will also be stable. 
Therefore, a coexistence equilibrium with two nearly identical prey has the same stability as the equilibrium 
for the corresponding one-prey systems. During damped oscillations onto a stable coexistence equilibrium, 
or diverging oscillations away from an unstable one, the clones will oscillate nearly in phase with each other 
and inherit the cycle period of the one-prey system. 



Top edge: The rightmost portion of the top edge also corresponds to the cusp in the coexistence equilibrium 
region, so the stability here is also the same as that of the (predator + vulnerable prey) system. In general, 
as \/k\ approaches the upper limit of the coexistence equilibrium region when p\> p* (the curved portion), 
the stability of the two-prey system approaches that of the (predator + defended prey) system with l/k\ 
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approaching 1 /k 2 . This must be stable if the (predator + vulnerable prey) system is stable, because the 
defended prey is always more stable, as noted above. If the (predator + vulnerable prey) system cycles, then 
there will be instability as pi — > 1 along the top edge. 

However, there is always stability near the top edge for p\ — > p*, as follows. Along the straight portion of the 
top edge (pi < p*), as 1 /k\ approaches the edge, the coexistence equilibrium converges to a limit with y = 0, 
while along the curved portion the limiting coexistence equilibrium has x 2 = 0. So near their intersection at 
pi = p*, both x 2 and y approach 0. Condition J20i then implies that the (predator + defended prey) system is 
stable, so the coexistence equilibrium is stable near the top edge for pi just above p*. By continuity, there is 
an open region of (p i,ki) values near p i — p*,ki = k 2 where the coexistence equilibrium is locally stable. If 
the (predator + vulnerable prey) system is only weakly unstable then this stability region may be quite large 
(Fig|5]panel B), but it cannot reach either the bottom or right edges. 



Left edge: Finally, consider the edge pi — 0. The steady states simplify to 



l-Z-x 2 -y, x 2 = Q 



, y = Q(m-l] 



(h-k 2 ) 

k\ + k 2 (m — 1) 



(29) 



where Z = ^-y . The coexistence equilibrium exists for # < ^- < where # is the value of 1 jk\ that solves 
x 2 + y + Z = 1 , noting that y depends on ki . The Jacobian matrix at d29t is 



J = 



— aixi —aixi — a\Xi 

-a 2 x 2 (-a 2 +gyF 2 )x 2 -(a2+gF)x 2 







where setting pi = and p2 = 1 gives F = j- 



Cli 



ghW 2 
and 

mkj 







(30) 



(ki+Z) 



(31) 



Near the lower limit of the left edge, we know that the system inherits the stability of the (predator + vulner- 
able prey) system. Above the lower limit we can use the Routh-Hurwitz criterion (Appendix|B) to determine 
stability. The coefficients cq and ci have common factor X2gF 2 > 0. Dividing this out gives modified coeffi- 
cients 



c ( ) = aixigk b F > 0, ci — k b (ci2 +gF) —aixi, C2 — a\x\ +X2(a2 ~gy~F 2 ) 



(32) 



and the stability conditions remain the same: cq, ci, C2 > 0, C1C2 > co. Extensive numerical evaluations of 
the coefficients as 8 is varied indicate that loss of stability occurs when the condition cic 2 — cq > is violated 
- the equilibrium is stable if this condition holds and unstable if it fails. Assuming this is true, loss of stability 
along the left edge occurs via a Hopf bifurcation (AppendixlBl. 
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5.4 The structure of evolutionary cycles 



The stability analysis above delimits the situations in which evolutionary cycles occur. As illustrated in 
Figure [5j\, they arise when the p\ versus k\ tradeoff curve passes (with decreasing p\) from the region 
of stable coexistence equilibria near p\ = p*,k\ = ki to the region of unstable coexistence equilibria with 
p\ i=s 0. For \jk\ below the dash-dot curve in Figure |5j\ the defended prey cannot invade the vulnerable 
prey-predator limit cycle (see Figure^. As 1 jk\ increases, the defended prey becomes persistent and then 
increases in average abundance. As \/k\ — * 1 jki the characteristic features of evolutionary cycles emerge: 
longer cycle period and out-of-phase oscillations in predator and total prey abundance. 

The phase relations on evolutionary cycles can be seen in the dominant eigenvector of the Jacobian matrix 
for the unstable fixed point (Figure^. There is a codominant pair of complex conjugate eigenvalues, and ( 
because det(7) < ) the third eigenvalue is real and negative. As the defended prey has very low palatability, 
the predator and the vulnerable prey have the classical quarter-phase lag. Here the phase angle is 90"; because 
eigenvectors are only defined up to arbitrary scalar multiples, including arbitrary rotations in the complex 
plane from multiplication by e' e , only the relative phases of eigenvector components are meaningful. As 
1 /k\ increases, the eigenvector components for the two prey types become out of phase with each other 
(« 180° phase angle, right column of Figure^}. As a result, the predator and total prey densities are out of 
phase with each other. 

In the next section we show that these phase relations become exact as the limit k\ — > ki is approached, for 
a general version of the model in which we do not specify the functional forms of the predator and prey 
functional responses. 

6 Evolutionary cycles in a general two-prey model 

In this section we analyze the limiting phase relations in evolutionary cycles, for p\ <C 1 and low cost to 
defense, without specifying the functional forms of the prey and predator functional responses. We consider 
a two-prey, one-predator model that (after rescaling) can be written in the form 

x i =x i (f(X,y,8 i )-p i yg(Q)),i=l,2 

y = y{Qg{Q)-d) (33) 

where as usual X — x\ +X2 is the total density of prey and Q = p \x\ + P2X2 is the total prey quality as perceived 
by the predator. The key assumption in (I33i is total niche overlap in the prey types (e.g., because they are two 
clones within a single species), which is reflected in / being a function of X. To model the trophic relations, 
/ is assumed to be strictly decreasing in X and nonincreasing in y, and h(Q) = Qg(Q) is strictly increasing in 
Q. Model (1331 includes the Abrams-Matsuda model (a two-prey version of the Rosenzweig-MacArthur 
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model with Lotka-Volterra competition between the prey) and the two-prey, one-predator chemostat model 
analyzed in the previous section of this paper. 

The parameter 0; in J33l > represents the prey-specific cost of defense, with / decreasing in 0. Because / is 
decreasing in X, we can parameterize / such that is the steady-state density for a single prey-type in the 
absence of predators, i.e. 

/(fl J 0,fl)=0. (34) 
As usual we number the prey so that pi < p 2 and therefore 9\ < 92, and rescale the model so that p2 = 1. 
Evolutionary cycles can be viewed as arising when 9\ f 02 with /'l < 1, such that when 9\ « 02 

1. There is a positive coexistence equilibrium with X\,X2 converging to positive limits while y — > as 

2. The coexistence equilibrium is an unstable spiral for 9\ w 02- 

At the end of this section, we show that Condition 1 always holds for 1331 . and in Appendix|F|we show that 
the coexistence equilibrium is always a spiral. To determine the limiting phase relations, we need to find the 
eigenvector corresponding to the dominant eigenvalue with positive imaginary part (see Appendix lAli: the 
relative phase angles of this eigenvector's components (in the complex plane) correspond to the phase lags 
between the corresponding state variables in solutions to the linearized system near the steady state. 

The Jacobian for J33I in the limit described above is 



h = 



xjx xxfx x\f y -p\x\g 

x%fx x 2 fx X 2 fy~X 2 g 





(35) 



The characteristic polynomial of J35i factors to show that the eigenvalues of J35i are fx(x~i +x 2 ) < and 
as a repeated root; the first has corresponding eigenvector (x\ ,JC2,0), and for there is the unique eigenvector 
(1,-1,0). The zero eigenvalue therefore has algebraic multiplicity 2 and geometric multiplicity 1. 

The long period of evolutionary cycles is explained by the fact that the dominant eigenvalues converge on a 
double-zero root as 9\ — > 02- The cycle period near the fixed point is inversely proportional to the imaginary 
part of the dominant eigenvalues, which converges to in this limit. 

To determine the limiting phase relations consider a small perturbation of the defended prey parameters, 
01 = 02 — £. For £ small, we show in Appendix|F]that the double-zero eigenvalue is perturbed to a complex 
conjugate pair of eigenvalues. To study cycles, we assume that these have positive real part. That is, near 
the double-zero root the (scaled) characteristic polynomial is perturbed to leading order from p{z) = z 2 to 
p(z) — (z — £a) 2 + eb 2 for some b > 0, so the perturbed eigenvalues have 0(e) real part and imaginary parts 
zky/ebi to leading order (here i = We need to determine the corresponding perturbed eigenvectors. Let 
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vvo denote the unperturbed eigenvector (1,-1,0), and let wq + w e be a perturbed eigenvector corresponding 
to the complex eigenvalue with positive imaginary part, scaled so that its first component is 1. The first 
component of w e is therefore 0. The perturbed Jacobian is Jq + eJ\ for some matrix J\ . Then 



(Jo + £Ji)(w + w e ) = y/ebi(w + w e ) + O(e). 
Using Jq\vq = and keeping only leading-order terms, gives 

Jow e = y/ebiwQ. 

Let w e = (0,W2,W3); then writing out (I37> in full, wi and W3 satisfy 



xi fx xify-pmg 

Xlfx X 2 fy-X 2 g 



Vi>2 
U-'3 



(36) 



(37) 



(38) 



H'2 



W) 



and W3 must be pure imaginary, because the unique solution to the real part of (I38i is (0,0). Writing 
(^/ebi)zj and solving for the z's, we find that zi < and Z3 > 0; specifically 



Z2 




(X\ +X 2 )fy - (j>lXl +X2)g 




'Xfy-d 


. Z 3 . 




-(xi+x 2 )fx 




_ -Xfx _ 



using the fact that (from the second line of j33» 

Qg(Q) = d. 

So to leading order the eigenvector corresponding to eigenvalue ^/ebi + o(^fe) is 



(39) 



(40) 



-\-^/eBi 



for some B,C > 0. 



(41) 



Now add total prey as a fourth state variable to the system. The corresponding eigenvector component is the 
sum of the first two components in J41i (see AppendixIXl: 



(42) 



1 

-\-y/eBi 
y/sCi 

We can multiply each component of J41> by an arbitrary real constant without affecting the phase angles, so 
we can consider instead (1,-1 — \/eBi, i, —i \ . Then as e — > the vector giving the relative phases for prey 
1, prey 2, predator, and total prey becomes 



[1 



1 i 



(43) 



The components of the limiting phase-angle vector J43i lie exactly on the coordinate axes. The two prey 
types (first and second eigenvector components) are exactly out of phase; the predator and total prey (third 
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and fourth components) are exactly out of phase; and there is a quarter-period lag between the vulnerable 
prey and the predator. This holds in the limit Q\ — > 02 for p\ < p* such that the coexistence equilibrium 
remains is an unstable spiral when 0\ < 62 ■ 

Finally, we now show that for pi sufficiently small there always exists coexistence equilibria as 0i f 62 such 
that xi,X2 approach finite limits while y — ► 0. Evolutionary cycles then occur whenever these equilibria are 
unstable. 

We must have 62 > Q - this is just the condition that the vulnerable prey, at steady state, is a sufficient supply 
of food for the predator to increase when rare. Fix a value of y > - we will show that for any y sufficiently 
small, there will be a value of 9i « 02 such that there is a coexistence equilibrium with y as the equilibrium 
prey density. A coexistence equilibrium must satisfy and the conditions 

f(X, y ,9 l )= Pi yg(Q). (44) 

Let X be defined as the solution of 

f(X,y,e2)=yg(Q). 

We have /(02,O, 62) =0, soX = 02 + 0(y). Then condition (1441 for prey type 1 can be satisfied by choosing 
01 = 2 + O(y) such that 

f(X,y,e 1 )=p l yg(Q) 
The last thing needing to be shown is that there exist x\ > Q,X2 > such that 

Xl+X2 = X, P\X\+X2 = Q. 

This is a system of two linear equations in the two unknowns X\,X2 whose solution (by inverting a 2 x 2 
matrix) is 



X\ 


1 


X-Q 


. *2 




. Q-Pi* . 



(45) 



For y small, X ss 2 > Q so *i > in J45I . And for pi sufficiently small, p\X <Q so x 2 > 0. So there do exist 
positive prey steady-states giving a coexistence equlibrium for y near 0, and their limiting value (as y — > 
and the corresponding 9\ — > 62) is given by ( 145 l l with X = (h- 

The argument above also identifies when p\ is "sufficiently small" for the limiting values of Jci,jC2 to be 
positive as Q\ — » 62- The limiting value of x\ is always positive, and the limiting value of x 2 is positive if 
p\X — Q > in the limit, i.e. if pi 62 < Q- Thus, the coexistence equilibrium region extends up to the line 
01 = 62 so long as pi < p* = (5/02- 
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7 Discussion 



The model studied in this paper is three dimensional, with a few fairly tame nonlinearities - just like the 
Lorenz equations. So it is not surprising that a complete mathematical analysis of the model has not been 
possible. Nonetheless, we have come a long way towards our goal of characterizing how and when rapid 
evolution can affect the ecological dynamics resulting from predator-prey interactions. 

Our primary questions concern the generality of the phenomenon of "evolutionary" limit cycles in predator- 
prey interactions, and the conditions in which such cycles might be observed. A combination of analysis and 
numerical studies suggests that evolutionary dynamics are not omnipresent, but neither are they knife-edge 
phenomena existing only in a narrow range of parameter values. Instead, the types of cycles observed by 
Yoshida et al. 14 II 1431 are both robust and general. They occur in a specific but substantial and biologically 
relevant region of the parameter space, and in a general class of predator-two prey models that includes a 
two-prey model with Lotka-Volterra prey competition terms [2 3"! 1251 . and the standard two prey chemostat 
model II 1II23I I41 1 with mechanistic modeling of resource competition between the prey. 

We have shown that evolutionary cycles arise through a bifurcation from a stable coexistence equilibrium, that 
occurs when defense against predation for the defended type remains relatively inexpensive but nevertheless 
becomes very effective. Cryptic population dynamics, where the predator cycles but the total prey density 
remains nearly constant, occur as a limiting case when effective defense comes at almost zero cost |43 1. These 
regions in parameter space are biologically relevant because empirical studies have shown that defense - be it 
against predation or against antimicrobial compounds - can arise quickly and can be both highly effective and 
very cheap [4 16 42 1. For example, Gagneux et al. 1 16 1 showed that in laboratory cultures of Mycobacterium 
tuberculosis (TB) mutants, prolonged treatment with antibiotics results in multi-drug resistant strains of TB 
with no fitness costs for resistance, and furthermore that most naturally circulating resistant TB strains are 
either low or no cost types. Indeed, fitness tradeoffs in the production of defensive structures and compounds 
are notoriously difficult to demonstrate, and in many empirical studies, no fitness tradeoff was actually found 
0|9]|37l. 

We close by listing some open questions. "Proving things is hard" (Hal Smith, personal communication), but 
others may succeed where we have not. Concerning the model in this paper, 

• When does the Jacobian at a coexistence equilibrium have a pair of complex conjugate eigenvalues? 
There will be 3 real, negative eigenvalues if the two prey types are very similar and the interior equi- 
librium for the (predator + vulnerable prey) exists and is a stable node. However, our numerical results 
suggest the full system (at a coexistence equilibrium) has complex conjugate eigenvalues whenever the 
(predator+vulnerable prey) system has an interior equilibrium with complex conjugate eigenvalues. 

• Can there be coexistence of the predator and both prey on a limit cycle or other attractor, even when 
there is no coexistence equilibrium? Numerical evidence suggests that the answer is "no" for the 
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chemostat model: for k\ below (above) the range of values at which a coexistence equilibrium exists, 
the defended (vulnerable) prey type outcompetes the other. As it is difficult to distinguish between per- 
sistence and slow competitive exclusion numerically, it is likewise hard to map reliably the parameter 
region where both prey coexist on a nonpoint attractor. 

• On the bifurcation curve < p\ < p* , k\ = kj, the Jacobian of the general model J33l > has zero as a 
double root with algebraic multiplicity 2 and geometric multiplicity 1 . Generically, this situation gives 
a Takens-Bogdanov bifurcation |26 1. Do the higher order conditions for Takens-Bogdanov, which hold 
generically, hold for our model {0? 

• A general two-prey, one-predator chemostat can exhibit a wider range of dynamic behaviors than we 
have observed in a system where the prey differ only in their p and k values (see [40 1 and references 
therein). Indeed, these predicted dynamics have been observed in other experimental systems |8|. 
The absence of some dynamics from our system, if true and verifiable, could indicate a qualitative 
difference between within-species evolutionary dynamics resulting from prey genetic diversity, and 
food-web dynamics with one predator feeding on a several prey species whose within-species heritable 
variation is much smaller than the functional differences among prey species. 

Finally, how robust are the phenomena of evolutionary and cryptic predator-prey cycles in more complex 
food webs involving multiple predator and prey species? 
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Appendices 



Appendices |A] and [S] summarize some general results useful to us here, and contain nothing original. In 
Appendix |C] we derive the expressions for coexistence steady states in the reduced and rescaled two-prey 
chemostat model, and in AppendixlDlwe derive the Jacobian matrix and prove that it has negative determinant 
at any coexistence steady state. In AppendixlElwe derive the conditions in which a limit cycle of the (predator 
+ edible prey) subsystem can can be invaded by the defended prey. Finally, in Appendix|F]we show generally 
that for realized cost Q\ sufficiently close to 62 and < p\ < p*, the coexistence equilibrium for the general 
model J33I always has a pair of complex conjugate eigenvalues. 



A Appendix: Eigenvectors and phase relations 

The contents of this Appendix appear to be well-known, but we have not seen them summarized anywhere in 
print. We consider oscillations in a linear system 

i = Jx (46) 
resulting from the real matrix J having complex conjugate eigenvalues 

A , A = a±ib with b > 0, 

where i — and the overbar denotes complex conjugation. The corresponding eigenvectors are also a 
complex conjugate pair w,w. The resulting oscillatory terms in solutions of J46> are of the general form 
Ae^'w + Be^'w. In order for these to be real (as solutions of ( 146 \ must be), we must have B =A. Then writing 
A = re' e , r > 0, the solutions are proportional to 

z(t) ee e ie e ibt w + e- ie e - ib 'w. (47) 

We are interested in the relative phases of the oscillations by different components in z(t ). Write wj = rie'^ 
for the f h component of w. The /* component of z(t) is then 

r ,( e K<t>j+e+bt) +e -i{pj+e+bt)) = 2r/cos(fy + + bt). (48) 

The relative phases of the /* and k th components in solutions proportional to z{t) is therefore given by <f>j — ^ . 
When this is near components j and k are oscillating in phase, and when it is near ±7T they are oscillating 
nearly out of phase. 

We are interested in the phase difference between the predators and total prey density. For that we can use a 
linear change of variables 
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In transformed coordinates the Jacobian matrix becomes AJA , and Jacobian eigenvectors w are trans- 
formed to Aw. The dominant eigenvector component for X\ + X2 is therefore the sum of the components for 
x\ and^2. 



B Appendix: Stability conditions 

In this Appendix we review criteria for local stability of equilibria in a three-dimensional system of ordinary 
differential equations. 

The diagonal expansion ([35 1, section 4.6) is an expression for det(A+D) where A is square and D is diagonal. 
For D =xl and A of order n it states that 

det(A +xl) = x" +x"- l Ti (A) +x"- 2 T 2 (A) + ■■■ + T„(A) (49) 

where 7) (A) is the sum of all principal minors of order j (a principal minor of order j is the determinant of 
a /' x j submatrix of A whose diagonal is a subset of the diagonal of A - that is, a submatrix obtained by 
selecting n — j diagonal elements of A and deleting the row and column containing that element). Note that 
T n (A) = det(A) and T x (A) = trace(A). 

For a 3 x 3 matrix the characteristic polynomial is 

p{X) = det(A/-A) = X 3 +c 2 X 2 + ciX+c . (50) 

Comparing with J49I and noting that and that Trj{— A) = {— 1) ; 7> ; '(A), we have 

Co = r 3 (-A) = -det(A), c l =T 2 (-A) = T 2 (A), c 2 = trace(-A) = -trace(A). (51) 

In the notation of i5Q\ . the Routh-Hurwitz stability criteria for order-3 systems (May 1974) is 

co > 0, c\ > 0, c 2 > 0, c\C2>cq. (52) 

Loss of stability through a Hopf bifurcation occurs when the third condition in d52t is violated, with the c, all 
positive 1 18 1. 



C Coexistence steady states for the rescaled chemostat model 

We consider here the two-prey model (II Oi Setting y = and solving gives the steady state value of Q, 
Q = -\. We solve 
x\ = X2 = imply 



= —K, We solve forX and y as follows. Defining Z = 1 — X — y and noting that , 8 - = i, the conditions 

^ g-l J b k h +Q Q> 



mZ Pl y mZ^_l = l 



h+Z Q k 2 +Z Q 
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Solving fl53b for y gives two expressions which remain equal within the coexistence region: 



Pi 



(m-\)Z-k\ 



ki+Z 

Setting the two expressions for y equal, we can solve for Z 

1 



(m-\)Z-k 2 



k 2 +Z 



2(l-p 1 )(m-l) 



C + VC 2 +4(m-l)(l- /?1 )2fe 1 fe 2 



(54) 



(55) 



where 



C = ki{l+pi(m-l))-k 2 ({m-l)+pi). 



Finally, recalling that Z — 1 — X — y, then X = 1 — Z — y. Expressions for x\ and x 2 in terms of X and Q are 
derived and shown in the text. 



D Jacobian at a coexistence equilibrium 



The general expression ill 2b for Jacobian entries at a coexistence equilibrium implies that all entries in the i th 
row of the Jacobian have common factor JEj, so det(7) = jcijC2ydet(/) where J(i,j) = with xi, = y. Let F 
denote the steady state per-capita feeding rate for the predator, 

1 



h+P\X\+p2X2 

and the a, are defined by d3 1 1 with Z = 1 — x\— X2~y; equation (1551 gives the general expression for Z. 
Taking the necessary partial derivatives, 



(56) 



J 



-ai+gp\yF 2 -a l +gp 1 p 2 yF 2 -a x -gpiF 
-a 2 + gp\p 2 yF 2 -a 2 + gpjyF 2 -a 2 - gp 2 F 
PighF 2 Pigk b F 2 



(57) 



We now show that the determinant of the Jacobian is always negative for the general model J33L and therefore 
for the chemostat model, unless p\ = p 2 . For d33l with the scaling p 2 = 1 we have 



fx - p\yg' fx -piyg' fy-pig 
fx-piyg' fx-yg' fy-g 
h'p h' 



(58) 



where g ~ g(Q),g ! — g'{Q) and n — h'(Q).h{Q) = Qg(Q). Then using basic products of determinants, det(7) 
equals 



h! 



fx fx fy-plg 
fx fx fy-pig 

P\ 1 



■ h! 



fx fx fy-Pll 

(pi-1), 
PI 1 



[1- Pl ) 2 h'gf x 



(59) 



which is negative (unless p\ = 1) because n > 0,g > and fx < 0. 
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E Appendix: Invasion of an edible prey limit cycle 



Following 1 3 1 we give here the condition for invasion of a predator + edible prey limit cycle by a rare defended 
prey type. Along the limit cycle we have (logy) = and therefore ^ kl+x 2 ) = 1- By Jensen's inequality, this 

& \_, l /m(l-x 2 -y) N (6()) 



implies that k g ty > > 1, and therefore (x 2 ) > Q- We also have (logjC2) = along the limit cycle, so 



^kb+x2/ \k 2 + 1 —xz — y 

A rare defended prey can invade if (logjq) > 0, i.e. if 

° < ( fl rT y l - *rf - - i ) = ( r ( i r T y l ) (&) - 1 

+ 1 -jc 2 -;y ^+x 2 / \fci + i-*2-)7 \h+x2/ 

Using J60l > and simplifying, we get the invasion condition in terms of p\,k\)\ 

^(CMTI' (61) 

where 

K; + 1 - X2 - y 

Note that the right-hand side of ( I6U can be computed for all using one long simulation of the (predator + 
vulnerable prey) system, and yields pi as a function of k\. 



F Appendix: Eigenvalues for Q\ ] Q21P1 < 



We show here that for Q\ sufficiently close to 9 2 and < p\ < p* in the general model 03l >. the coexistence 
equilibrium always has a pair of complex conjugate eigenvalues. As Q\ — > 02. m this range of /?i values 
y — > 0, so we set y = e <C 1 and use a series expansion in e of the characteristic polynomial (i.e. we regard Q\ 
as a function of y with all else held fixed, rather than vice versa). The Jacobian at the coexistence equilibrium 
is an 0(e) perturbation of (I35> and so to leading order has the form 



J(e) 



A + ea u A + £a\2 B + Efln 

C + £fl21 C + £«22 D + EC123 
£fl31 £fl32 



(62) 



with A,B,C,D < 0, and 031 = 771032 > (the last holds because y/y is a function of Q = p\x\ +X2 with the 
scaling p 2 = 1). J(0) has eigenvalues zero (with algebraic multiplicity 2) and A +C < 0, and we need to 
approximate the near-zero eigenvalues for e small. The characteristic polynomial of J(e) is a cubic in X but 
the near-zero eigenvalues are at most 0(\fe), so for our purpose the A 3 terms in the characteristic polynomial 
can be neglected. This leaves a quadratic polynomial in X, which will have complex conjugate roots if its 
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discriminant is negative. Using Maple to compute the characteristic polynomial of J62t . discard A 3 terms and 
expand the remainder about e = 0, to leading order in e the discriminant is 

4e (a 32 - a 3 1 ) (A + C) {AD - BC) 

which will be negative if AD — BC > 0. Referring to J35I some algebra gives 

AD-BC = x l x 2 f x g{ Pi -l) 

which is positive because fx < 0, as desired. 
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A: Single clone dynamics 



B: Two clone dynamics 



a> 

CD 

> 
B 

"D 



LO 
O 




CD 
_3 

CO 
> 

CO 
-I — I 

CD 
T3 



LO 

o 




' , WA'AWA'AWA'AWA'AWA'AWAWA'A'A'A'w'a'a'a'a 

^A'A'A'A-A'A'A'A'A'A-A'A'A-A'A'A'A'A'A'A'A-A-A'A'a'aW 
a' a'a' a *a 'a' a' a'a' a'a - 



AWAWA'A'A'A'A'A'AWA' 
'A'A'A'A'A'A'A'A'AWA'A'A' 
'A'A'AWA'AWA'A'A'A'A'A' 



trait value, p 



i 1 1 1 1 r 

0.0 0.2 0.4 0.6 0.8 1.0 
trait value, p 



C: Single clone, no age structure 



D: Two clones, no age structure 
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Figure 2: Dynamics of one- and two-prey models as a function of the palatability p of the defended prey type and the 
dilution rate S. Panels A and B show results for the "full" model including predator age-structure; panels C and D show 
results for the "reduced" model without predator age structure. A,C. Dynamics of a single prey system. Open circles 
show predator-prey cycles, and white space indicates equilibrium. B,D. Dynamics of a two-prey system. The model is 
scaled so that the undefended prey type has p=l. Filled circles indicate that both types coexist and cycle together. Open 
circles show short predator-prey cycles with only the defended type (p < 1). Cross-hatching indicates the defended and 
undefended prey coexisting at a stable equilibrium. Open triangles indicate equilibrium between predator and defended 
prey and white space indicates equilibrium between predator and vulnerable type. In the model with age structure (panel 
B), predator extinction occurs for S > 1.5 
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One-prey simplified model 
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Figure 3: Bifurcation diagram for the rescaled, reduced clonal model with one prey type. 



A: Edible prey cycles B: Edible prey stable 
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Figure 4: Prey coexistence equilibria. The shaded gray regions indicate parameters (pi,l/k\) for prey type 1 giving 
a coexistence equilibrium (stable or unstable) with the vulnerable prey type p2 = 1. At S = 1.5 (panel A) prey type 2 
cycles, and at 8 = 2.0 (panel B), prey type 2 is stable. The dashed lines show a representative tradeoff curve J4j, assuming 
roughly 50% reduction in growth rate as the cost of being 100% defended. Here k c = 0.017, a\ = 1.0, and «2 = 0.0165. 
In panel A the dash-dotted line indicates where the defended type can invade the limit cycle of the predator and vulnerable 
prey type. 
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A: Vulnerable prey cycles 
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B: Vulnerable prey barely cycles 
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C: Vulnerable prey weakly stable 
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D: Vulnerable prey stable 




Figure 5: Stability of coexistence equilibria for the reduced two-prey model. In each panel, the horizontal axis is the 
palatability p\ of the defended prey with the model scaled so that p2 = 1. To remain consistent with Abrams (cf. |3J, 
Figures 1-3) the vertical axis is 1/ki, scaled so that and 1 correspond to the lower and upper limits of the coexistence 
equilibrium region (Figure[4}. The Jacobian matrix and its eigenvalues were evaluated at an even 50 x 100 grid of values. 
Lighter gray indicates that the equilibrium is stable, darker gray that it is unstable; in all cases the computed eigenvalues 
with largest real part are a complex conjugate pair. The vertical black line is at p\ = p* , the value where the straight and 
curved segments of the upper limit of the coexistence equilibrium region meet. The dashed curve in panels A and B is the 
tradeoff curve k\ = k c + 0:2(1 — Pi) a ' , with k c = £2 = 0.054, (Xi = 1; (X2 = 0.05 at p\ values lying within the coexistence 
equilibrium region. The dash-dotted line represents the minimum 1 jk\ values at which the defended prey can invade 
the (predator + vulnerable prey) limit cycle (see Appendix|EJ. Parameter values for these plots are as follows: panel A., 
S — 1.5; B., S = 1.75; C, S = 2.0 and D., 5 = 2.1. All other parameters are unchanged and are as shown in Table 2. 
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Fraction of defended prey: Average (defended/total) 
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Figure 6: Contour plot of the long-term average fraction of defended prey. The horizontal axis is the palatability p\ of 
the defended prey, with the model scaled so that p2 = 1 . The vertical axis represents 1 /k\ , with and 1 corresponding to 
the lower and upper limits of the coexistence equilibrium region (Figure[4}. Numerical solutions of the model were used 
to compute the long-term average value of x\/(x\ +X2) f° r parameter values such that the (predator + vulnerable prey) 
system (same parameter values as panel A of Fig. [5}. In the lighter-gray region the coexistence equilibrium is stable. In 
the darker-gray region the equilibrium is unstable. The vertical black line is at p\ = p* , the value where the straight and 
curved segments of the upper limit of the coexistence equilibrium region meet. Parameter values are as in Table 2 with 
8 = 1.5. 
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Figure 7: Coexistence of edible and defended prey on a limit cycle. Parameter values for all plots were 5 = 0.9, m = 
3.3/5, g = 2.3/5, k 2 = 0.05, k h = 0.2, p 2 = l,p\ = 0.08. Values of k\ were 0.4 (top row), 0.1 (center row) and 0.055 
(bottom row). In each row the leftmost panel shows the dynamics of total prey and predator densities, the center panel 
shows the dynamics of the two prey types, and the rightmos^anel shows the phases of the Jacobian matrix components: 
l=def ended prey, 2=edible prey, 3=predator, 4=total prey. 



